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Abstract — Thermal underwater glider can harvest thermal 
energy from the temperature difference between surface sea 
water and deep sea water. PCM (phase change material) which 
is filled in underwater glider’s heat exchange tube solidifies and 
melts periodically with volume contraction and expansion when 
the thermal underwater glider ascends and descends between 
the surface and deep sea water, thus realizing energy conversion 
and storage. In this paper the numerical method of equivalent 
specific heat capacity was used to stimulate the phase change 
process. Using the calculated transient temperature field, we can 
get the position of solid-liquid moving boundary and mass 
fraction of liquidus PCM versus time. The theoretical time for 
complete melting was much less than the experimental time, it 
was because of the existence of residual air within the 
phase-change energy storage system. The thermal conductivity 
of liquidus PCM needs to be modified, the correct coefficient 
was 0.7-0.9. Combining with the nonlinear model of 
phase-change energy storage system, system pressure versus 
time was draw. The curve showed good consistency with 
experimental data. Related parameters were identified, air 
solubility in liquidus PCM was 0.05, and residual air within the 
system was 0.01-0.03.The residual air had great influence on the 
melting time of PCM and stored energy of phase-change energy 
storage system. The correct coefficient of thermal conductivity 
was one key factor to determine the curve of system pressure 
versus time and its value was fixed by comparing theoretical and 
experimental melting time of PCM. The fixed value of thermal 
conductivity was empirical and needs to be further studied. 

Index Terms —Underwater glider, Phase change, Residual air, 
Thermal energy 


I. INTRODUCTION 

Phase-change energy storage system is a prospective 
technique for the storage of thermal energy by using phase 
change material (PCM). During its phase change, PCM can 
absorb or release large amount of latent heat. This technique 
is wildly used in domestic and industrial domain like 
buildings, green houses, solar energy storage and etc. 

PCM has been a main research topic for the recent 20 years 
[l]-[3]. In 2004, P. Lamberg et al. invested the melting and 
freezing processes in phase change material storage 
numerically and experimentally [4]. In 2005, H. Kumano et al. 
studied direct contact melting with hydrocarbon mixtures as 
the PCM [5]. In 2012, H. H. Al-Kayiem et al studied thermal 
behavior of encapsulated phase change material energy 
storage [6]. Mostly, PCM is used as thermal energy storage 
medium and its volume change is neglected or avoided. The 
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total volume of PCM contracts when it changes from liquid to 
solid phase and expands reversely. This feature can be 
utilized in phase-change driving mechanism like actuator and 
micro valve. L. Klintberg et al developed a large stroke, high 
force actuator which used paraffin as phase change material 
[7]. R. Liu designed single-use, thermally actuated paraffin 
valves for microfluidic applications [8]. D. C. Webb et al 
designed the first thermal underwater glider using 
phase-change driving mechanism [9]. Q. Kong et al studied 
the phase change process for underwater glider propelled by 
ocean thermal energy numerically and experimentally 
[10]-[11]. Z. Hong-Wei et al putted forward the application 
and improvement of the interlayer thermal engine powered by 
ocean thermal energy in an underwater glider [12]. 

In this paper, combining the numerical method of 
equivalent specific heat capacity and previously established 
nonlinear model of energy storage process, the air solubility 
in liquidus PCM and residual air within the system was 
identified according to experimental data. The theoretical 
curve of system pressure versus time was gained and showed 
good accordance with experimental curve. The experimental 
data also verified the numerical method indirectly. 

II. Procedure for Paper Submission 

The whole configuration of thermal underwater glider is 
shown in Fig. 1, the heat exchange tube which is filled with 
PCM is attached to the main body, and PCM can absorbs and 
releases heat through the tube wall with outside sea water 
during phase change process. The cross section of heat 
exchange tube is described in Fig. 2, the middle of heat 
exchange tube is rubber hose which is filled with oil and 
connected with hydraulic circuit. The remaining space of heat 
exchange tube is filled with PCM of constant mass. When 
PCM solidifies, oil can be absorbed into the rubber 
diaphragm from the oil tank with low pressure; when PCM 
solidifies, oil can be extruded into the accumulator thus 
increasing the system pressure and realizing energy 
conversion and storage. 

In order to analysis the system pressure versus time during 
melting process, it needs to combine equivalent specific heat 
capacity method and the previously established nonlinear 
model of energy storage process for thermal underwater 
glider. Through the numerical method, we can get mass 
fraction of liquidus PCM versus time, and through the 
nonlinear model of energy storage for thermal underwater 
glider, we can get the system pressure versus mass fraction of 
liquidus PCM. After that, we can draw the curve of system 
pressure versus time and analysis the deviation between 
theoretical and experimental data and do the following 
parameter identification. 
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Fig. 1 Schematic of thermal underwater glider 


A. Equivalent specific heat capacity method 

During the solidification process of PCM, the outer layer of 
PCM solidifies continuously, and the boundary of solid and 
liquid phase moves toward the center point of the cross 
section continuously. During the melting process, the 
condition is similar. This problem is called moving boundary 
problem or Stefan problem [13]. Except for some simple 
situation which can be solved by analytical method, in actual 
engineering application, this problem can only be solved 
using numerical method [14]-[15]. 

In this paper, the melting process is simulated using the 
numerical method of equivalent heat capacity. Based on the 
calculated transient temperature field, the mass fraction of 
liquidus PCM versus time is gained consequently. 

The melting process of PCM in heat exchange tube is a 
three-dimensional heat transfer process, and the solving zone 
also changes because of the density difference of liquidus and 
solidus PCM. In the following numerical simulation, the 
solving zone is assumed to be constant and meshed using 
fixed grid not moving grid for simplicity [16]. In order to 
simplify the numerical method, assumption is made as 
follows: 

1) Ignore the heat resistance of heat exchange tube; 

2) Axial temperature distribution of heat exchange tube is 
equal; 

3) Phase change material is isotropic and homogeneous; 

4) Ignore the effect of temperature on density and other 
thermal parameter of PCM; 

5) Influence of system pressure on phase change process is 
also neglected. 


Thermal Energy 


Thermal Energy 


O 


otl ~^§! 


os 


liquid PCM 

Fig. 2 Cross section of heat exchange tube 

Based on the above assumption, the original problem can 
be simplified as one-dimensional heat transfer problem. 
Under cylindrical coordinate system, the energy control 
equation can be described as 

57 . ,d 2 T 1 c>7\ 


pc — = k (—— h-) 

dt dr 2 r dr 


(1) 


For a control volume with constant mass, assume the phase 
change happens in a narrow temperature range 
[T m -AT,T m +AT] 

capacity is given by 


then the equivalent specific heat 


C +C, 


c = - 
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( 2 ) 


For a control volume, mass fraction of solidus phase is used 
to describe different zone and is defined as follows 
1 T <T m - AT 


fs 


T-T 


2A T 
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T>T m + AT 
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Where T s =T m — AT . Its value range is [0 1]. When f s = 0 , 
it means the control volume is solid phase, when f s = 1, it 
means the control volume is liquid phase, and when 
0 < f s < 1, it means mash zone with solid and liquid phase. 
For a control volume of constant mass, mass fraction f s is 
also the volume fraction of solid phase. 

Since the convection of liquidus PCM is remarkable during 
melting process [17]-[18], its effect should be taken 
consideration by introducing effective thermal conductivity. 

= (). 18 * Rp 5 (4) 

Using this numerical method, the transient temperature 
field can be calculated. Through the following process, the 
curve of mass fraction of liquidus PCM versus time in the heat 
exchange tube can be got. 


1. Tracking the moving boundary of solid and liquid 
phase 

For a meshed solving zone, when the temperature of two 
adjacent nodes meets the constraint 
T(i)>T m >T(i + \) 

boundary by 


,then we can get the position of moving 


: r(i+1) + E r(l+1) 


(5) 


7 ( 1 )- 7 ( 1 + 1 ) 

Where r(i+l) , r(i) are the position of nodes I,i+1 in radial 


coordinate system respectively and r m the position of moving 
boundary. 


2. Mass fraction of liquidus PCM during melting 
process 

Since the remaining space of heat exchange tube which 
is filled with PCM is annulus, and the melting process is 
simplified as a one-dimensional heat transfer problem, mass 
fraction of liquidus PCM during phase change process can be 
calculated by dividing the liquidus area with the whole 
annulus area. 


w = 



( 6 ) 


Where R0 is the outer radius of the annulus area, and Ri the 
inner radius of the annulus area. 


3. Modification of thermal conductivity of liquidus PCM 
When comparing the theoretical time for complete melting 
with experimental time, the former is much less than the latter, 
the main reason is residual air within the system. It forms 
bubble in the beginning of melting process thus decreasing the 
melting process severely. The actual thermal conductivity of 
liquidus PCM is less than the theoretical one. It is needed to 
be modified by correction coefficient beta. The formula of 
beta is given by 
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T 

beta = — (7) 

T 

With T t the whole time for complete melting in theory, and 
T e the whole time in experiment. 


B. Nonlinear model of phase-change energy storage system 


Thermal underwater glider can harvest environmental 
energy from the temperature difference of surface sea water 
and deep sea water. The actual component includes heat 
exchange tube, connected with hydraulic circuit, the inner oil 
tank, accumulator, etc. It can be simplified as a phase-change 
energy storage system which is showed in Fig. 3, 4. The 
performance of the phase-change energy storage system is 
significant to harvest and conversion of thermal energy from 
the sea water. Because the system is inevitably mixed with a 
certain amount of air, and even a small amount air will also 
have a severe effect on the process of increasing pressure 
phase. Based on the previous theoretical analysis, the formula 
for solving system pressure at different stage is obtained. 

The critical pressure when residual air within the system 
completely dissolves is described by 


P c = -^-*P 0 (8) 

s 0 * w 

During increasing pressure phase 1, the system pressure is 


-B + ^B 2 - 4A * C 
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During increasing pressure phase 2, the system pressure is 
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At the turning point of increasing pressure phase 1 and 
increasing pressure phase 2, we can get the constraint 
equation Pl=P2=Pc, namely 

s + a 
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Where 
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III. EXPERIMENT AND ANALYSIS 

Hexadecane is selected as PCM for thermal underwater 
glider because of its suitable melting temperature and certain 
amount of density difference. Other features of hexadecane 
are as follows: it is compatible with the heat exchange tube, its 
performance is stable for long duration and it is safe without 
toxicity or corrosivity. Detail property of hexadecane is listed 
in Table 1. 


Table 1. Property of PCM [5] 


Property 

Unit 

Value 

Melting temperature T m 

°C 

16-18 a 

Latent heat of fusion L 

KJ/kg 

236 

Thermal conductivity k s /ki 

W/mK 

0.313/0.14 

Specific heat c s /ci 

J/kgK 

1.64/2.09 

Density p J P/ 

Kg/m 3 

864/773 


a The actual melting temperature of PCM is a constant value, it 
is transformed into a narrow range which is required by the 
numerical method. 

The research team of Tianjin University had done a series 
of experiments on Qiandao Lake of Zhejiang province in 
September 2013.The temperature distribution of the testing 
water is shown in Fig. 5. The surface water temperature is 
about 27 °C, and the bottom water temperature is about 10.5 
°C. At the beginning of testing, underwater glider was set 
aside on the bottom of the lake for a long enough time. After 
PCM in the heat exchange tube completely solidifies, rapidly 
raises it to the surface and set aside. PCM begins to melts with 
volume expansion. Since oil was extruded from the rubber 
hose of heat exchange tube to the accumulator, the system 
pressure increases continuously until PCM completely melts. 
Underwater glider automatically records the system pressure 
versus time during the whole melting process. The possible 
melting and temperature change of PCM during underwater 
glider ascending from bottom to surface water is ignored since 
time is little. 


1) When the value of residual air is zero, air solubility in 
liquidus PCM is 


[ l-a-—) + 1 K- 1 

^0 PsO 


B*T, 


(15) 


P *T 

r i 1 o 


2) When the value of air solubility is already calculated, we 
can get the value of residual air by 


A. Raw experimental data 

The curve of system pressure versus time can be divided 
into three sections roughly as shown in Fig. 6 and Fig. 7. 
Taking Fig. 7 as an example, after thermal underwater glider 
ascends to the surface water, PCM begins melting. Because 
the existence of residual air within the system, system 
pressure increases slowly until it reaches the value of critical 
pressure, and this stage is named as increasing pressure phase 
1. When system pressure exceeds the critical pressure, 
residual air dissolves completely in melted PCM, and system 


36 


www.erpublication.org 


















Mathematical and experimental analysis of energy storage for thermal underwater glider 


pressure increases sharply until it reaches the set value of high 
pressure tank, this stage is named as increasing pressure phase 
2. After that, system pressure keeps constant as the set value 
of high pressure tank and oil is transferred into high pressure 
tank continuously until PCM melts completely, this stage is 
named as energy storage phase. These three stages are divided 
by two turning points. 



Fig. 5 Temperature distribution of testing water 



Fig. 6 System pressure versus time for group #1~#6 



Fig. 7 System pressure versus time for group #7 


In this paper, only increasing pressure phase 1 and 
increasing pressure phase 2 were analyzed, and mainly focus 
on the turning point of those two stages. When simulating the 
melting process using the numerical method of equivalent 
specific heat capacity, the influence of system pressure on 
phase change is ignored, namely the property of PCM keeps 
constant no matter what the value of system pressure is. 

The experimental data is showed in Table 2. Among the 8 
groups, the first 7 groups are raw experiment data, group #8 is 
calculated by mean of corresponding item respectively. 

B. Related parameter identification 

Taking group #8 as an example, using the experimental and 
theoretical value (Tt=2909) of whole time for PCM 
completely melting, then correct coefficient beta can be 
calculated. After modifying the thermal conductivity of 
liquidus PCM, the numerical method of equivalent specific 
heat capacity was applied again to simulate the phase change 
process. The curve of mass fraction of liquidus PCM versus 
time was drawn. 

In experimental curve, after finding the turning point and 
its corresponding time, we can get the corresponding mass 
fraction of liquidus PCM in the above simulated curve. It 
means although PCM melts continuously, the system pressure 
didn’t show obvious increase until certain amount of PCM 
melts because of the existence of residual air in the system. 

According to the above two steps, calculate the 
corresponding mass fraction of liquidus PCM for group 
#1~#7, as shown in Fig. 8.When assuming the typical value of 
air solubility is 0.08, calculated results showed residual air for 
all groups are negative, indicating the actual value of air 
solubility is less than 0.08. Here, according to the calculated 
value of mass fraction in #1~#3 team, assume the value of 
residual air for these three team is zero, then by Equation (15), 
we can get the value of air solubility for these three teams. The 
calculated values show good consistency. For the following 
#4~#8 team, using the already calculated value of air 
solubility and Equation (16), we can get the value of residual 
air respectively. The detail results of air solubility and 
residual air for all the teams are shown in Table 3. 



Testing group 

Fig. 8 Correct coefficient and mass fraction at turning point 
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Table 2. Raw experimental data 


Testing group 

1 

2 

3 

4 

5 

6 1 

8 a 


413 

334 

348 

347 

371 

358 341 

3594 

Time for complete melting (s) 

5 

1 

6 

2 

9 

5 9 







851 

997 952 

697. 

Time for turning point (s) 

527 

435 

457 

664 



6 






12 

12 12.4 

12.1 

Steady system pressure (MPa) 

12 

12.1 

12.2 

12.5 



7 

a Group #8 is not actual raw experimental data, the corresponding item is 

calculated by average. 


Table 3. Parameter identification. 








Testing group 1 2 

3 


4 

5 

6 

7 

8 

0.052 0.052 0.052 

0.052 

0.052 

0.052 



air solubility s 0 2 3 

5 


3 

3 

3 

0.0523 

0.0523 




0.010 

0.016 

0.027 



residual air a 0 a 0 a 

0 a 


7 

5 

7 

0.0277 

0.0111 


a Group #1~#3 are assumed to have zero value of air solubility. 


C. Theoretical curve of system pressure versus time 

Taking group #8 as an example, based on the correction 
coefficient beta and numerical method, we can draw the curve 
of mass fraction of liquidus PCM versus time during PCM 
melting process. In the beginning PCM melts quickly, then its 
melting rate gradually slows down. The melting rate 
decreases continuously during the melting process. The whole 
time for complete melting is 3557s using modified thermal 
conductivity of liquidus PCM. 

Based on the two identified parameters of air solubility and 
residual air within the system, the curve of system pressure 
versus mass fraction of liquidus PCM is obtained. When the 
mass fraction of liquidus PCM is 0.482, it changes from 
increasing pressure phase 1 to increasing pressure phase 2, 
and the system pressure increases sharply after the turning 
point. 

Combining these two curves, we can draw the curve of 
system pressure versus time. As shown in Fig 9, when the time 
is 697.2s, it changes from increasing pressure phase 1 to 
increasing pressure phase 2, and the system pressure at this 
turning point is 0.254MPa. The turning point and changing 
trend of this curve showed a very good consistency with 
#1~#7 experiment data. It also verified the numerical method 
and nonlinear model indirectly. 



IV. Conclusion 

The mass fraction of liquidus PCM versus time during 
melting process is gained using the numerical method. 
Combining the nonlinear model of phase-change energy 
storage system, the curve of system pressure versus time is 
obtained and showed a good consistency with experiment. It 
also verified the numerical method indirectly. The identified 
parameters are as follow: 

1) During the melting process of PCM, residual air within 
the system can severely decreases the melting rate, resulting 
the whole time for complete melting in experiment much 
more than that in theory. It needs to modify the thermal 
conductivity of liquidus PCM. The value range of correction 
coefficient beta in this paper is 0.7-0.9; 

2) The step to identify the relevant parameters is 
established. The calculated value of air solubility in liquidus 
PCM is 0.05, and the value range of residual air within the 
system is 0.01-0.03; 

The correct coefficient of thermal conductivity is empirical 
and needs to be further studied. 

NOMENCLATURE 

p Density (kg/m 3 ) 

c Specific heat (J/kg • K) 

L Latent heat of fusion (J/kg) 

k Thermal conductivity (W/ m • K) 

/ Paritcle fraction of liquid/solidified PCM 

T Temperature (K) 

t Time(s) 

r Radial coordinate (m) 

L Latent heat (kJ/ kg) 

T m Melting temperature (K) 

R Rayleigh number 
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Diameter of heat exchange tube (m) 
Diameter of rubber diaphram (m) 
Correction coefficent 
Mass fraction of liquid/solidified PCM 
Air solubility under standard condition 
Pressure (Pa) 

Equivalent 
Liquid phase 
Solid phase 

refference / atmosphere 
Low 
High 

Phase change material 
Phase-change energy storage 
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